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ABSTRACT 

An analytical approach to the problem of replenishing task forces 
at sea is investigated, using a random walk model. The type of 
replenishment operations considered consists of two replenishment 
ships of different types, each replenishing L combatant ships of 
different types. The combatant ships are initially distributed into 
queues of length M and N between the two replenishment ships. As each 
combatant ship completes its first replenishment, it then enters the 
queue of the other replenishment ship. Thus each combatant ship is 
replenished by both replenishment ships. Replenishment times for the 
combatant ships are assumed to be independent, positive valued random 
variables whose distributions are dependent upon both replenishment 
ship type and combatant ship type. The replenishment of the task 
force is completed whenever 2L replenishments, two for each combatant 
ship, are completed. The total replenishment time is shown to be the 
maximum of certain partial sums of the individual replenishment times. 

A general expression for the distribution function of the total re- 
plenishment time is sought by various methods of analysis . A computer 
simulation is used to gain a better understanding of the form of the 
distribution function of total replenishment time for a specific under- 
way replenishment operation. Some examples are discussed. 
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I. THE UNDERWAY REPLENISHMENT OPERATION 



The underway replenishment operation is the means by which needed 
supplies are provided to Navy fleet units while underway at sea. It 
provides these units with the capability of operating at great distances 
from Navy shore facilities, enables them to remain on station for ex- 
tended periods, extends their period of effective combat readiness, 
and limits diversion from their assigned mission. The supplies pro- 
vided may include, but need not be limited to, petroleum products, 
missiles, ammunitions, stores, and provisions. However, no matter 
what is to be transferred between the transferring replenishment ships 
and the receiving combatant ships, several points are apparent, the 
underway replenishment operation increases the vulnerability of the 
engaged units, restricts their mobility, and diverts the involved 
combatant ships from their principal missions. The Navy is continually 
striving for the development of techniques and methods which will 
reduce the limitations imposed on fleet operations by the necessity 
for underway replenishment. Principally, the Navy desires to reduce 
the vulnerability of the engaged units, increase the unit mobility 
while engaged, reduce the time required to complete replenishment, 
and limit interference with the assigned missions of the combatant 
units . 

To attain these goals many improvements of a material nature have 
been accomplished. The introduction of faster replenishment ships, 
such as the combat stores ship (AFS), reduces the vulnerability and 
increases the mobility of the engaged units by allowing underway 
replenishments to be conducted at higher speeds than before. The 
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introduction of multi -commodity replenishment ships, such as the fast 
combat support ship (AOE) and the fleet replenishment tanker (AOR) , 
reduces the total number of replenishments required since several needs 
can be filled during a single replenishment. The introduction of the 
helicopter as a transfer agent affords greater separation of task force 
units, reduces the demands for alongside replenishments, and minimizes 
interference with the assigned mission of the involved combatant ship. 

The introduction of new wire-highline transfer systems, such as the 
fast automatic schuttle transfer cargo handling system (FAST), provides 
a higher rate of cargo transfer and an increased load capacity, thereby 
reducing the required time alongside. In addition, internal cargo 
handling improvements in both the receiving and transferring ships 
have resulted in reduced breakout and strikedown times, thus again 
resulting in reduced replenishment times. All these improvements have 
had one element in common. Their accomplishment has resulted in a 
reduction in the actual replenishment time for an individual ship, and 
a subsequent reduction in the total replenishment time for the involved 
task force. 

Improvements have not been limited to those of a material nature. 

The underway replenishment operation has also been viewed as a problem 
in operations analysis. The objective of such an approach is a better 
understanding of the underway replenishment operation as viewed in an 
analytical manner. 

The underway replenishment operation can be viewed in many analytical 
contexts, but more readily lends itself to representation as a multi- 
server queueing model with a finite number of customers. The number of 
servers can be established by the number of replenishment ships involved. 
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while the number of customers can be established by the number of ships 
in the task force scheduled for underway replenishment operations. It 
should be noted that, under normal conditions, a definitive sequence of 
individual replenishments has been established prior to the commencement 
of the underway replenishment operation. This is done as a method of 
optimizing replenishment ship utilization, while hopefully minimizing 
total replenishment time for the task force. Thus, while the underway 
replenishment operation may be viewed in terms of a queueing problem, 
it has definite job-shop scheduling problem overtones. 

One aspect of a typical queueing model must be ignored for the 
underway replenishment operation. It is normal to consider all service 
times to have come from the same probability distribution in a typical 
queueing model. Since many different types of replenishment and com- 
batant ships are involved in a typical underway replenishment operation, 
this seems to be an unreasonable assumption for any realistic model of 
the operation. Analysis of actual underway replenishment data has 
supported this contention and has yielded approximations for the 
appropriate distribution functions for various typical underway 
replenishment operations . 
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1 1 . THE UNDERWAY REPLENIS HMENT MODE L 



The basic underway replenishment model that was developed by 
Waggoner [1] and later expanded by Patterson [2] is the domain of 
inquiry for this thesis. This model assumes that two replenishment 
ships and L combatant ships are scheduled to conduct underway replenish- 
ment operations; a total of 2L replenishments are to be conducted, L by 
each replenishment ship. Prior to the commencement of the replenishment 
operation, the L combatant ships are given a specific replenishment 
sequence and assigned to queues of length M and N, where M + N = L, 
associated with the first and second replenishment ships, respectively. 
For each combatant ship this means that a definitive order of replenish- 
ment has been established and that each combatant ship has initially 
been assigned to one of the two queues. In addition, this means that 
the individual combatant ship both succeeds and precedes the same ships 
in each queue and that this queueing sequence is inalterable. 

The movement of the combatant ships during the underway replenishment 
operation is envisioned as a sequence of steps in a random walk. This 
random walk is performed on an L x L grid on which the coordinates of 
each grid node represent a possible combination of the numbers of com- 
batant ships replenished by each replenishment ship, i.e., the node 
(i,j) on the grid indicates that i combatants have been replenished by 
the first replenishment ship and j by the second replenishment ship. 

It is obvious from the replenishment sequence that certain grid nodes 
represent infeasible combinations; it is infeasible, for example, for 
the first replenishment ship to have completed M + k replenishments 
unless the second replenishment ship has completed at least k, where 
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k = This infeasibility can be stated in constraint equations. 

If the abscissa, Z^, of the grid denotes the number of replenishments 
completed by the first replenishment ship, while the ordinate, Z 2 , of the 
grid denotes the number of replenishments completed by the second re- 
plenishment ship. The following equations constitute the necessary 
constraints : 



Z, = 



Z, = 



z, = 



Zr. = 



Zr. = 



0 

M + Z. 
L 
0 

N + Z, 



Z^ = L 



All feasible nodes either lie on or within the boundaries established 
by the intersection of these constraint equations. Figure 1 depicts 
an L X L grid with the constraint equations superimposed. 

In order to establish a common terminology for later discussion 
and investigation of the model, various definitions and theorems 
stated by Patterson [2], are presented. 

Definition 1 

The CONSTRAINT BOUNDARY of the random walk grid is defined as follows 
(see Figure 1): 

a. the UPPER BOUNDARY consists of those nodes along the line 
Z 2 = L with abscissa k, where M < k < L-1; 

b. the RIGHT BOUNDARY consists of those nodes along the line 
Z^ = L with ordinate h, where N < h < L- 1 ; 

c. the LOWER BOUNDARY consists of those nodes along the line 
Z^ = M + Z 2 , where 0 < Z 2 < N-1; and 
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Underway Replenishment Model 
Random Walk Grid 
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d. the LEFT BOUNDARY consists of those nodes along the line 
Z 2 = N + Z^, where 0 < < M-1. 

Definition 2 

A STEP is a horizontal or vertical line between adjacent nodes of the 
grid; a PATH is a sequence of non-negative steps which starts at the 
origin and terminates at (L,L) and does not cross the constraint 
boundary; and and denote the times required to complete the ith 
horizontal and jth vertical step, respectively, of any path. 

From Definition 2 it is evident that is the ith replenishment 
time for the first replenishment ship, while Y^ is the jth replenish- 
ment time for the second replenishment ship. 

Definition 3 

The set of all feasible paths from the origin to (L,L) is decomposed 
into four mutually exclusive and exhaustive subsets: 

A:, the set of all paths that pass through the node (L-1,L) and 
do not pass through any nodes on the lower boundary; 

B: the set of all paths that pass through the node (L-1,L) and 

pass through at least one node on the lower boundary; 

C: the set of all paths that pass through the node (L,L-1) and 

do not pass through any nodes on the left boundary; and 

D: the set of all paths that pass through the node (L,L-1) and 

pass through at least one node on the left boundary. 

In addition, subsets B and D can be further decomposed: 

Bj^: the subset of all paths in B which last pass through the 

lower boundary at node (M + k,k) for k = 0,1,..., N-1 

D^ : the subset of all paths in D which last pass through the 

left boundary at node (h , N+h) for h = 0,1,..., M-1. 
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Theorem 1 



The total replenishment time is the actual time required to complete a 
path through the grid and is equal to one of the following values: 

a. the total length of time required for any path in A is 

L 

T = E X. ; 

A i=l 1 

b. the total length of time required for any path in is 

L k+1 

T = E + E Y. for k = 0,1,. . . ,N-1; 
i=Hfk+l j = l ^ 

c. the total length of time required for any path in C is 

L 

T = E Y ; or 

^ j=i " 

d . the total length of time required for any path in is 

, h+1 L 

= E X. + E Y. for h = 0,1,..., M-1. 
i=l j=N+h+l 

Waggoner and Patterson both assumed that the and Y^ were inde- 
pendent, random variables from the family of exponential probability 
distributions. Waggoner tacitly assumed that the replenishment rate 
was dependent only upon the replenishment ship type and that replenish- 
ment rates X and \Jl could be used for all replenishments conducted by 
the first replenishment ship and the second replenishment ship, res- 
pectively. Patterson generalized this assumption to the extent that 
the replenishment rate was allowed to be dependent upon both the 
replenishment ship and the combatant ship types involved. Thus re- 
plenishment rates X^ and were assumed for the ith replenishment by 
the first replenishment ship and the jth replenishment by the second 
replenishment ship, respectively. It should be noted here, however, 
that Patterson developed Theorem 1 without referral to his assumption 
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of exponentiality . Thus Theorem 1 holds without any assumption being 
made regarding the distribution of the individual replenishment times. 
An interesting comparison can be made regarding the investigative 
approaches of Waggoner and Patterson, Waggoner approached the model 
by considering the first step taken and the subsequent sequence of 
steps of shortest duration taken during the random walk. Thus the 
exponentiality assumption, by means of the memoriless property, 
insured that the duration of each step continued to be exponentially 
distributed. This greatly simplified his derivation of the Laplace 
transform for the distribution function of the total replenishment 
time. On the other hand, Patterson reversed the sequence and con- 
sidered the last step taken and the subsequent sequence of steps of 
longest duration in the reverse order. In this case, as noted above, 
the exponentiality assumption was not applied, Patterson, then using 
the exponentiality assumption, was able to develop the Laplace trans- 
form for the distribution function of total replenishment time only 
for the case of two combatant ships initially distributed one to each 
replenishment ship . 
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III. TOTAL REPLENISHMENT TIME 



The previous section has intimated that total replenishment time 

can be one of certain feasible values. It was, in fact, stated that 

total replenishment time was either T . , or T_ (for some integer k 

A D 

where 0 < k < N-1), or T , or T ^ (for some integer h where 0 < h < M-1). 
In this section, by means of the following theorem and proof, it is 
shown that total replenishment time is actually the maximum of those 
f eas ib le va lues . 

Theorem 2 

The total replenishment time, i.e., the time required to move from node 

(0,0) to node (L,L) of the random walk grid, is the maximum of all 

feasible values. Let T denote total replenishment time, then 

T = max (T. ; for k = 0,1, •• • ,N-1; T_ ; for h = 0,1, .. . ,M-1). 
Ad L D 

Proof 

a. Assume that the path taken was from A. Thus it must be shown 
that for k = 0,1,..., N-1; ; and for 

h = 0,1, . . . ,M-1. 

1) Because the path is in A, it must not touch the lower 
boundary. Thus it must reach the line Z 2 = k + 1 (see Figure 1) before 
it reaches the line = M + k for any k = 0,1,..., N-1. This implies 
that 

M+k k+1 

E X. > E Y. for any k = 0,1,...,N-1. 
i=i ^ j=i ^ 

L 

Now by adding the partial sum E X. to both sides of the above 

i=M+k+l ^ 

inequality, it is seen that T > T for any k = 0,1,..., N-1. 

A D 
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2) Because the path is in A, it must reach the upper boundary 
before it reaches the right boundary. This implies that 



3) Let the path in A first touch the line = N + h (see 
Figure 1) at the node (r,N+h) where h < r < L-1 for any h = 0,1,..., M-1. 
Since the line Z^ = N + h must be reached vertically, the (r + l)th 
horizontal step must have begun before the node (r,lH-h) is reached. 

The (IH-hd-l)th vertical step is either just beginning (h+1 < r < L-1) 
at node (r,N+h) or will begin on completion of the (r + l)th horizontal 
step ( r = h). Since the upper boundary must be reached before the 
right boundary, this implies that either 



L 



L 




Therefore, ^ • 



i=r+l ^ 



L 

S X. > 



j=N+h+l 



L 

S Y. if h+1 < r < L-1, or 

U-V.J-1 J 



L 

E X. > 



L 

E Y. if r = h. 

u.Ka.1 J 



i=r+2 ^ 



j=N+h+l 



In either case, the above inequalities imply that 



L 



L 



E X. > E Y. for any h = 0,1,..., M-1. 
i=h+2 " j=lH-h+l J 



h+1 



Now if the partial sum E is added to both sides of the last 



inequality, it can be seen that T >T for any h = 0,1,..., M-1. 

iJ 



b. Assume that the actual path taken was from Bj^ where k is fixed 

k k n 

and 0 < k < N- 1 . Thus it must be shown that T_ > T. ; T_ > T for • 

BABB 



n k, 0 < n < N-1; > T ; and > T ^ for h = 0,1, .. . ,M-1. 

“ ” B U B D 



1) Because the path is in Bj^, it must reach the line = M + k 



(see Figure 1) before reaching the Z^ = k + 1. Thus 
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k+1 >H-k 

S Y. > S X. . 
j=l J i=l ^ 

L 

Nov? by adding the partial sum Z X. to both sides of the above 

i=lH-kd-l ^ 

Ic 

inequality, it can be seen that > T^. 

2) Because the path is in it must pass through the node 

(lH-k,k). Thus the (Mfkfl)th horizontal step begins at the same time 
as the (krh2)th vertical step. Since the path cannot pass through any 
node (M+u,u) where k < u < N-1, it must reach the line Z 2 = u + 1 
(see Figure 1) before reaching the line = Mfu. These conditions 
imply that 

M+u ^ u+ 1 

E ^i ^ Z Y. , for k < u < N-1. 
i=M+k+l j=k+2 J 

L k+1 

Now by adding the two partial sums E X. and H Y. to both sides 

i=lH-u+l ^ j=l J 

k u 

of the above inequality, it can be seen that T_ > T_ for k < u < N-1. 

B B 

3) Let the path in Bj^ first touch the line Z 2 = v+1 (see 
Figure 1) at the node (s,v+l) where 0 < s < M+v and 0 < v < k. The 
(v+2)th vertical step begins at this node, while the (s+l)th horizontal 
step either just begins (s = Mfv) or has already begun (0 < s < Mfv-1). 
Since it is necessary that the line Z 2 = krhl (see Figure 1) not be 
reached until after the line Z^ = M + k is reached, this implies that 
either 

k+1 m-k 

EY.> EX. ifs=M+vor 

j=v+2 ^ i=s+l ^ 

k+1 lifk 

EY.> EX. ifO<s<M + v-1. 

j=ry+2 ^ i=S+2 ^ 

In either case, this implies that 
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k 



kfl Mfk 

S Y.> S X. forO<v<k. 
j=v+2 i=Mfv+l ^ 

v+1 L 

Now if the partial sums S Y. and S X, are added to both sides of 

j=l ^ i=lH-k+l ^ 

k V 

the last inequality, it can be seen that T > T for 0 < v < k. 

U — 

4) Because the path is in the (k+2)th vertical step and 
the (M+k+l)th horizontal step begin at the same time. Since the 
upper boundary must be reached before the right boundary, this implies 
that 



S X. > i: Y.. 

k=l«-k+l ^ j=k+2 J 

kfl 

Now by adding the partial sum S Y. to both sides of the above 

j=i J 

inequality, it can be seen that T_ > T _ . 

5) Because the path is in it cannot reach the line = k+1 

(see Figure 1) before reaching the line Z^ — Mfk, which implies that 



k+1 



IH-k 



S Y. > S X. , 
■ i=l " 



j=l J 



This in turn implies that 

k+1 h+1 

(1) S Y. > E X. for any h = 0,1, . . . ,M-1. 
j=l J i=l ^ 

Further, let the path in Bj^ first reach the line Z 2 = N + h (see 
Figure 1) for any h = 0,1,..., M-1 at some node (w,N+h) where M+k < w < L. 
Since the path reaches = N + h through a vertical step, the (N+h+l)th 
step is just beginning at node (w,N+h) while the (w+l)th horizontal step 
has already begun. Since the upper boundary must be reached before the 
right boundary, it is necessary that 



L L 

E X.> E Y. forM+k<w<L. 
i=w+l ^ j=N+h+l ^ 



21 



I 

I 

I 



9 \ 



This in turn implies that 

L L 

(2) S X. > S Y. for any h = 0,1,. . . ,M-1. 
i=Mfk+l " j=W-h+l J 

Now if the inequalities (1) and (2) above are added together, it can be 
k h 

seen that T > T for h = 0,1,..., M-1, 

D D 

c. If the path taken was in C, is larger than the other 
feasible values. This is proven analogously to part a. 

d. If the path taken was in , for some h = 0,1,..., M-1, 

is larger than the other feasible values. This is proven analogously 
to part b. 

This concludes the proof of Theorem 2. 

By way of illustrating the result of Theorem 2, an example where 

constant individual replenishment times are involved is discussed. 

It is assumed that = A and = B for i,j = 1,2,...,L for the 

general case of L combatant ships, distributed into queues of length 

M and N, associated with the first and second replenishment ships, 

respectively. If neither N = 0 nor M = 0, the following are the 

feasible values for T: 

T = LA 
A 

= (N-k)A + (k+l)B for k = 0,1, . . . ,N-1 

D 

Tc = 13 

= (h+l)A + (M-h)B for h = 0,1, . . . ,M-1. 

T, the maximum of all feasible values, is clearly equal to either LA 
or LB, depending upon the values of A and B. If A > B, then T = LA ; . 
if B > A, then T = LB. If N = 0, all feasible paths are in D and the 
following are the feasible values of T: 

Td*’ = (h+l)A + (L-h)B for h = 0,1,..., L-1. 
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T is clearly equal to either (A + LB) or (lA + B), depending upon the 
values of A and B. If A > B, then T = (LA + B) ; if B > A, then 
T = (A + LB). If M = 0, the preceding argument holds if the set D 
is replaced by the set B. 
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IV. IDENTIFICATION OF THE GENERAL PROBLEM 



In the preceding section it has been shown that the actual total 
replenishment time, denoted T, is the maximum of all feasible values. 

Thus, while the exact composition of T as expressed in terms of the 
individual replenishment times is known, the general form of its dis- 
tribution function as expressed in terms of those for the individual 
replenishment times remains unidentified. In regard to this last 
statement, in this or any succeeding section the following notation 

has been used: the distribution function of X. is denoted F. , while 

1 1 

that of the Y. is denoted G. ; if densities are assumed, they will be 
J J 

denoted by the appropriate small letter with subscript. 

The composition of T and the distribution function of T when L< 2, 
while of little practical interest, are presented below as an introduction 
to the analytical techniques used in the attempt to identify the general 
distribution function of T, when L > 3: 

a. if L = 1 with M = 1 and N = 0, then T = and its dis- 

tribution function is the convolution of F^ and G^; 

b. if L = 2 with M = 2 and N = 0, then T = X^ + max (X^jY^) + Y^ 
and its distribution function is the convolution of F^, ^ 2 ^^, and G 2 , 
where H 2 ^ denotes the distribution function of max (X 2 >Y^), i.e., 

F 2 G 1 ; 

c. if L =2 with M = 1 and N = 1, then T = max (X^,Y^) + max (X 2 ,Y 2 ) 
and its distribution function is the convolution of and ^ 22 ^ where 
H. . denotes the distribution function of max (X.,Y.), i.e., F.G. for 

i = 1,2. 
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It should be noted that for all cases investigated, either above 

or beyond this point, an assumption that M > N has been made. The 

symmetrical nature of the model has insured that any results obtained 

here are valid upon the interchange of M with N, with Y ^ , and 

with G. for the case when M < N. 

J 

The general distribution function of T can be represented by the 
following probability statements; 
a. if N = 0, then 



^+1 

P(T < t) = P Z + 



L 

S ^ 
j=hH-l 



< t 



for h 




b . if N > 1 , then 

, L L k+1 

P(T<t)=p(zX.<t; Z X.+ Z <tfork = 0,l,...,N-l; 
i=l ^ i=M4-k+l ^ j=l 

L h+1 L . 

Z Y. < t; Z X. + Z < t for h = 0,1, . . . 
j=l ^ i=l ^ j=N+h+l 

Viewing the summation terms involved in the two probability state- 
ments as random variables, denoted as T , it can be noted that when 

n 

N = 0 L such random variables are present, while if N > 1 L + 2 such 

random variables are present. Since it was desired to use the technique 

of transformation of variables to simplify these probability statements, 

it was obvious that additional random variables had to be introduced in 

order that the one for one transformation demanded by the technique 

could be attained. Since the transformation was to be from the X. and 

1 

the Y. into the T , L artificial random variables would have to be 
J n’ 

introduced if N = 0, while L - 2 would suffice if N > 1. These 
artificial random variables were selected to fulfill two purposes: 
primarily, they insure a one-to-one transformation, while secondarily. 
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they have no effect on the appropriate probability statements. For 

L 

the case of N = 0, random variables of the form Z X. for m = 1,2,...,L 

i=m ^ L 

are introduced. These summations are specifically less than Z X. + Y . 

L i=l ^ ^ 

For the case of N > 1 random variables of the form Z X. for 

L i=m 

m = 2,3,..., M and Z Y . for n = 2,3,..., N are introduced . These 

j=n ^ L L 

summations are specifically less than Z X. and Z Y. respectively. 

i=l ^ j=l ^ 

Thus for either case the introduction of these particular summations 
into the appropriate probability statement have no effect. 

The whole set of new random variables is defined below: 
a. if N = 0, 



T = Z X. for m = 1,2,. ..,L 
m . 1 

i=m 



h+l L 

. L/j 

1=1 j=h+l 



b. if N > 1, 



L 

T = E X. for m = 1,2,. . . ,M 
m . 1 

i=m 



^M+k+1 



L k+1 

E X. + E Y. for k = 0,1, ...N-1 
i=M+k+l ^ j=l ^ 



L 

Ttj. = S Y. for n = 1,2, . . . ,N 
U-n .1 > j j 

j=n 

h+l L 

1 = 2 X + E Y for h = 0,1,...,M-1. 

The Jacobian matrices of these transformations are shown in Figures 2, 
3 and 4. The value of their determinants is the subject of Lemma 1. 
Lemma 1 



The Jacobian determinant of the transformation matrices shown in Figures 
2, 3 and 4 is equal to 1, unless L = 2 with M = 1 and N = 1, i.e., if 
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JACOBIAN MATRIX OF T TRANSFORMATIONS WHEN N = 0 
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FIGURE 2 
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JACOBIAN MATRIX OF T TRANSFORMATIONS WHEN N = 1 
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FIGURE 3 
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JACOBIAN MATRIX OF T TRANSFORMATIONS WHEN N > 2 
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FIGURE 4 



|j| denotes the determinant of the Jacobian transformation matrix J, 
then 1 j| = 1. 



Proof 

a. Rows 1 through L + N have all 0*s to the left of the main 
diagonal elements and the main diagonal elements are equal to 1. 

b. Rows L + N + 1 through 2L-1 can be reduced to similar form by 
the application of the formula 

L+N+i L+N+i 1 1+1 ’ ^ ^ 



where R^ denotes the ith row before application of the formula, and 

R* . afterwards. 

1 

c. Row 2L can be reduced to a similar form by application of the 
formula 

'"’2L= ''2L - if N = 0, 

R* 2 l = ^2L " ^1 ' ^L+1 ^IH-2 ■ ^2 ^f ^ = 1, or 

*^*2L ^2L ■ ■ ^L+1 \+2 ^f N > 2. 

The case for L = 2 with M = 1 and N = 1 was specifically excluded 
from Lemma 1 since the determinant of its Jacobian transformation 
matrix equals zero. This results from the inherent dependence of the 
involved summations. It is noted that this is the only case in which 
no artificial random variables can be introduced. Of course, the 
distribution function of T in this case was given earlier in this 
section. 

A second transformation of variables was made at this point to gain 
further simplification. The following random variables were defined: 

V = T - T_^, for m = 1,2,...,2L-1 
m m m+1 

V = T 
2L 2L* 
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The determinant of the resulting Jacobian transformation matrix equals 
1 since the matrix is triangular with all main diagonal elements equal 
to 1, i.e.j all elements below the main diagonal are 0*s. 

At this point it was decided to assume that the distribution 
functions of the individual replenishment times were continuous and 
possessed non-vanishing density functions when the or are greater 
than 0. The joint density function of the and the can thus be 
written 






f (x^, . . . , . . ,y^) = “ 



TT f (x . ) n g.(y.) x^ > 0 and y^ > 0 



i=l ^ ^ j=l J J 



V 



for i J=l,2, . . . ,L 
otherwise , 



since it was previously assumed that the X^ and the Y^ were independent 
random variables. 

Since the appropriate Jacobian determinants are equal to 1, the 
joint density function of the can also be written in the same 



manner 



TT f. (x.) TT g.(y.) X > 0 and y. > 0 

ii=l j=l J J ^ J 

=1 for i, j = l,2, . . . ,L 

0 otherwise. 

The region in V-space for which g(v^^,V 2 , • . • ,V 2 j^) does not vanish can be 

represented in terms of the if the transformations relating the 

and the Y. to the V are used. These are given below: 
j m 

a. if N = 0, 



X. = V. for i = 1,2, . . . ,L-1 
11 

2L 

X = s V. 

L . T 1 
i=L 
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1 



b. if N = 1, 



c. if N > 2, 



Yj=V.^^+V^. for j = l,2,...,L-2 



2L 

r = y V + V 
L-1 . ^ i 2L-1 

i=L 



2L-1 

Y = - E V. 
^ i=l " 



X. = V. for i = 1,2, . . . ,M-1 
11 

2L 

X = - V. + + E V. 

L 1 L . - 1 

i=Li-2 



Y. = V. + . for j = l,2,...,L-2 

j j L+j 

’'l-I - ''l + ''m Vl + ''2L-1 

Yi^ = - 2V^ - V V. - + V 

1=2 



X. = V. for i = 1,2, . . . ,M-1 

11 J J > 

~ ^L+k+l for k = 0,1, . . . ,N-2 

\-l = ^ \-l + \+N 

2L 

X = - V, + + E V. 

^ ^ ^ i=L+W-l " 



Y. = for j = 1,2,... ,N-1 

J U-j 

^m-h=Vl + Wh for h = 0,l,...,M-2 

\-l = \ Vl + ^2L-1 
M 

Y=-VV-V +V 
L . , i L+1 2L 
1=1 



The above statements thus allow the original probability state- 
ments for the distribution function of T to be represented in terms 
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L 



I 






of the V in the following form: 
m 

2L 



P(T < t) = P E < t for m = 1,2, .. . ,2L^ 



i=m 



provided that the joint density function for the V does not vanish. 



m 



As an example, the joint density function of the V , expressed wholly 



m 



in terms of the V^, is presented below for the case N > 2: 



^M-1 



L-2 



. "/ i '“i \ Vi ^L- 1 

1=1 i=M 



2L N-1 



L-2 



fL<-v^+Vj+ E V^. ) n g,(v^_^.) TT 



i=Ixt-N+l 



j+l-N L+j- 



g(v^.V2,...,V2L> = 



M 



® L- 1 ^''m'*^L+i‘^'' 2L- 1 i L+ 1‘‘‘^2 

1 = 1 

if condition (1) below exists 
^0 otherwise. 

Condition (1) is expressed below: 

V . > 0 i = 1, . . . ,M-1 



''i + Vi+l > 0 i = 
''i + “n V» ^ “ 



2L 

-v^ + V + E V. > 0 

‘ i=u-m-i ‘ 



'' L+j > 0 J = 1 N-i 

''j+l-N + “L+j ^ “ j=N,...,L-2 

“m + “l+1 “2L-1 ° 

M 

- E V. - V + V > 0. 

. 1 1 L+1 2L 

1=1 



It is evident that the above discussion has not concluded in a 
truly desirable result, namely a general form of the distribution 
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function of T that would facilitate ease of computation. The resulting 
forms, those containing the and the do, however, exhibit prop- 

erties that may prove useful in future attempts in this area. As the 
matter now stands, not even a relatively simple problem, such as that 
when L = 3 with M = 2 and N = 1, can be directly resolved. Appendix A 
depicts this particular problem as an illustration of the general dif- 
ficulties encountered. In retrospect, it can be stated that the 
dependence re lationships between the T^, or those between the V^, 
impose restrictions that cannot be overcome without some convenient 
manner of handling the joint distribution function of dependent random 
variables, such as that developed in this section. 
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V. A SIMULATION OF THE UNDERWAY REPLENISHMENT OPERATION 



The preceding section, in which analytical methods were used, 
failed to produce the distribution function of T in a form feasible 
for specific computations. It was therefore decided to study the dis- 
tribution of T by computer generation of T realizations for a specific 
case of underway replenishment. The primary intent of such a study was 
to gain a better understanding of the form of the distribution function 
of T for typical underway replenishment operations. 

The general method used made direct use of the results of Theorem 
2, i.e., that T is the maximum of all feasible values. The underway 
replenishment operation was not directly simulated. Instead, the 
individual replenishment times were simulated and the resultant values 
of T calculated, using the results of Theorem 2. 

For the specific underway replenishment operation, it was decided 
to use an attack carrier task force replenishment operation. The task 
force was selected to be composed of eight combatant ships, an attack 
aircraft carrier (CVA) , a light guided missile cruiser (CLG), two 
guided missile frigates (DLG) , and four destroyers (DD) . The two ship 
replenishment group was selected to be composed of a fleet oiler (AO) 
and an ammunition ship (AE) . This particular selection was made because 
approximations for the distribution functions of the individual replenish- 
ment times for the above type ships were presented by Besecker in [3]. 

Thus Erlang distribution functions were assumed for the individual 
replenishment times; the specific gamma parameters, denoted k and \, 
are shown in Table I. The general form of the Erlang distribution 
function is given in Appendix B. 
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PARAMETERS FOR ERLANG DISTRIBUTION FUNCTIONS 

FOR 



INDIVIDUAL REPLENISHMENT TIMES 





REPLENISHMENT SHIPS 


AO 


AE 


k 




k 




c 














0 




CVA 


11 


0.09 


3 


0.03 


M 


s 












B 


H 


CLG 


7 


0.11 


2 


0.02 


A 


I 












T 


P 


DLG 


7 


0.11 


2 


0.02 


A 


S 












N 




DD 


4 


0.05 


3 


0.06 


T 















TABLE I 

In accordance with the model, the AO and AE were designated the 

first and second replenishment ships, respectively. The AO is then 

associated with the X., while the AE is associated with the Y.. It 

J 

was decided to establish a constant replenishment sequence at the AO. 
This sequence is CVA, DLG, DD, DD, CLG, DLG, DD, and DD. The replenish- 
ment sequence at the AE would depend upon the M and N combination being 
used for the generation of a T realization. The feasible sequences 
are shown in Table II. 

The individual replenishment times were simulated by the inverse 
transformation method discussed in Appendix B. The property that an 
Erlang random variable is the sum of k random variables, each expon- 
entially distributed with parameter A., facilitates this method. To 
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insure that the random variates generated were from the assigned dis- 
tributions, chi-square goodness of fit tests were made at the 0.01 
level of significance. This method is discussed in Appendix B. For 
all 16 distributions of the individual replenishment times, the null 
hypotheses that the generated random variates came from the appropriate 
distributions were of course accepted. 







REPLENISHMENT 


SEQUENCE AT 


AE 






M 


8 


7 
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3 


2 


1 


0 


N 
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1 
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3 


4 
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6 


7 
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CVA 


DD 


DD 


DLG 


CLG 


DD 


DD 


DLG 


CVA 


E 

P S 


DLG 


CVA 


DD 


DD 


DLG 


CLG 


DD 


DD 


DLG 


L E 
E Q 


DD 


DLG 


CVA 


DD 


DD 


DLG 


CLG 


DD 


DD 


N U 
I E 


DD 


DD 


DLG 


CVA 


DD 


DD 


DLG 


CLG 


DD 


S N 
H C 


CLG 


DD 


DD 


DLG 


CVA 


DD 


DD 


DLG 


CLG 


M E 
E 


DLG 


CLG 


DD 


DD 


DLG 


CVA 


DD 


DD 


DLG 


N 

T 


DD 


DLG 


CIjG 


DD 


DD 


DLG 


CVA 


DD 


DD 




DD 


DD 


DLG 


CLG 


DD 


DD 


DLG 


CVA 


DD 



TABLE II 

A special feature was provided in the implementation of the com- 
puter program to enable study of the dependence of T on the M and N 
combination being used. The program provided that upon the nth gen- 
eration of a T realization for a particular M and N combination, all 
ships would have the same individual replenishment times at both re- 
plenishment ships as for any other M and N combination, i.e., on the 
50th generation of a T realization for M = 3 and N = 5, the CLG would 
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have the same replenishment times at both the AO and AE as for the 50th 
generation of a T realization for any other M and N combination of the 
combatant ships. 

Five hundred T realizations were generated for each of the nine 
feasible M and N combinations. The random walk for a typical T 
realization is shown in Figure 5 as a random walk within the constraint 
boundary of the L x L grid and in Figure 6 as a random walk in two- 
dimensional- -time space. The T realization obtained by performing the 
random walk agrees with that obtained through an application of 
Theorem 2. 

For each M and N combination, the generated data were used to con- 
struct frequency histograms of T realizations. Figure 7 is a typical 
frequency histogram. From these histograms, one basic observation 
could be made: the frequency histograms showed definite positive 

skewness, i.e., they were skewed to the right. The obvious contention 
at this point was that the distribution function had the same character- 
istic shape as did those of the individual replenishment times, i.e., 
the shape of a gamma distribution function. 

A second observation could be made once the T sample means for the 
various M and N combinations were calculated. The results are similar 
to those cited by Milch and Waggoner in [4] for the expected values of 
T when the and are exponentially distributed with parameters \ 
and [i, respectively. The minimal T sample mean occurred when M and N 
were equal and a range of M and N combinations exists in which there 
is only a slight variation of the T sample mean values. The T sample 
means are plotted against M and N combinations in Figure 8. 

It was decided to investigate the contention that the distribution 
function of T was from the family of gamma distributions. Estimates of 
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FIGURE 5 
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necessary gamma distribution parameters k and denoted k and X, were 
obtained by the method of moments as shown in Appendix B. The basic 
null hypothesis to be tested would be that the T realizations for a 
particular M and N combination came from a gamma distribution with 

A A 

parameters equal to k and The chi-square method was choosen for 

testing goodness of fit. Despite the generally acceptable results, 
eight out of nine null hypotheses being acceptable at the 0.01 level 
of significance, it was evident that better results could be achieved. 

By studying superimposed frequency histograms of the observed T 
realizations and the expected T realizations under the gamma distri- 
bution of the null hypotheses for both the accepted and rejected cases, 
two observations could be made: 

a. the peak ordinate value for the number of observed T realizations 
was higher than that of the peak ordinate value for the number of 
expected T realizations for all M and N combinations; and 

b. the mode, the absissca for the peak ordinate value, for the 
number of observed T realizations either occurred to the right of or 
approximately coincided with the mode for the number of expected T 
realizations for all M and N combinations. 

A typical superimposed frequency histogram is shown in Figure 9. 

From these observations, the following generalization could be 
made : 

the true density function of T has a maximum value that occurs 
above and to the right of a gamma density function with parameters 

A A 

k and \ . 

Thus it was decided to use a method of estimating gamma parameters 
that utilized the maximum value of the gamma density function. This 
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method, termed the peak value method, is developed in Appendix B. 
Estimates of k and X, as determined by this method, are denoted k* and 
X*, respectively. Success was readily achieved by this method. The 
basic null hypothesis, that T for a particular M and N combination was 
distributed by a gamma distribution with parameters k* and X*, was 
acceptable at the 0.01 level of significance by the chi-square method 
for all nine M and N combinations. The results, including those 
obtained by the method of moments, are shown in Table III. It could 
be concluded that the distribution function of T could be successfully 
approximated by a gamma distribution function with parameters k* and X*. 



RESULTS OF METHOD OF MOMENTS AND PEAK VALUE 


ESTIMATION 


M/N 


T 


^t' 


A 

k 


A 

X 


k* 


X* 


8/0 # 


13. 12 


4.24 


40.58 


3.09 


46.19 


3.52 


7/1 


12.28 


4.08 


36.95 


3.01 


41.32 


3.36 


6/2 


11.83 


3.76 


37.24 


3.15 


44 . 63 


3.77 


5/3 


11.52 


3.18 


41.76 


3.63 


48.66 


4.23 


4/4 


11.49 


3.17 


41.67 


3.63 


47.70 


4.15 


3/5 


11.50 


3.19 


41.46 


3.60 


46 .34 


4.03 


2/6 


11.53 


3.22 


41.28 


3.58 


46.60 


4.04 


1/7 


11.95 


3.97 


35.91 


3.00 


45.06 


3.77 


0/8 


13.46 


4.71 


38.43 


2.86 


46.93 


3.49 



#: M and N combination that resulted in rejection of null 

hypothesis under method of moments testing 

TABLE III 
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APPENDIX A 



ILLUSTRATION OF PROBLEM WHEN L = 3 WITH M = 2 AND N = 1 



The elements of this particular problem are presented below in the 
context of Section IV: 

a. Considering the feasible paths, this problem involves random 
variables of the form: 



= Xi + X2 + X3 



4 = 



X 3 + Yi 



= 



Y1 + Y2 



T° = X 



+ Y. 



“ >‘1 ’‘2 



+ Y 
+ Y 
+ Y 



3 

3 

3' 



The distribution function of T can be written 

P(T < t) = P(X^+X2+X3 < t, X3+Y^ < t, Y^+Y2+Y3 < t, 

X^+Y2+Y3 < t, X^+X2+Y3 < t). 



b. After the addition of the appropriate artificial random variable, 
the problem involves random variables of the form: 



T, = 



X^ + X2 + X3 



= 



X 2 + X 3 



= 



X3 + Y1 



T, = 



Y + Y 
12 



T,. = X, 



+ Y. 



T, = 



X^ + X 2 



+ Y 
+ Y 
+ Y 



3 

3 

3* 
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The distribution function of T can be written 



P(T < t) = P(T. < t for i = 1,2,..., 6 ). 

c. After the transformation from the T to the V , the problem 

n m 

involves random variables of the form: 

= Xi 

V2 = X2 _ -Y^ 



^^3 = 



X, 



-Y -Y 
2 3 



V 4 =-x^ 



+Y, 



V3 = -X2 



+Y„ 



Ve = x^+X2 



+Y 3 . 



The distribution function of T can be written 



P(T<t) =P SV t for i = 1,2,. .., 6 ). 
m=i 

d. Considering the inverse transformations from the V to the X. 
and the , the problem involves the forms: 

Xi= 



X2 = + V2 



+ V, 



X 3 . -V 3 



+ V, 



+ Vc + V, 
5 6 



^1= 



+ V, 



Y = V + V 

I2 V2 



+ V, + V_ 
4 5 



Y^ =-2V - V 

3 1 2 



V, 






e. The joint density function of the V^, as wholly expressed in 



terms of the V , can be written 
m 
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'f 1 (V 1 ) f 2 (V ^+V2+V^ ) f 3 ( - V i+V3+V3+v^ ) x 



g(Vi,V2,V3,v^,V3.v^) = 



gl(Vi+v^)g2(v^+V2+v^+V3)g3(-2v^-V2-v^+v^) 

if condition (1) below holds 
otherwise 



Condition (1) is as follows 
V, 



V + V 
1 2 



V, 



V + V 
1 2 



> 0 

+ V, >0 



V + 
3 ^ 



+ V, 



+ V- + V, 
4 5 



V + V > V 
5 6 1 



> 0 
> 0 



V > 2V + V 
6 1 2 



+ V, 
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APPENDIX B 



RELEVANT THEORETICAL AND STATISTICAL FORMS 



1. The Gamma Distribution 



The gamma distribution function is given by 



F(x) = 



, k k"l -Xx , 

X X e dx 

r(k) 



for X > 0 and k > 0, 



/ k" 1 *"X 2 

X e dx. Let |j, m, and a denote the mean, the mode, 

o 

and the variance of the gamma distribution, then 

\l = k/X ; 

m = (k- 1)/X ; and 

2 i,a2 

a = k/X . 



If k is a positive integer, then the resulting gamma distribution is 
also termed an Erlang distribution. The Erlang distribution can be 
shown to be that for the sum of k identically distributed independent 
random variables from an exponential distribution with parameter X. 

If k is of the form n/2 with n a positive integer and X equals 1/2, 
then the resulting gamma distribution is also termed a chi-square 
d istr ibution . 

2 . Determination of Erlang Random Variables by the Inverse Trans- 
formation Method 

Let F(x) denote the distribution function for a random variable X 
and let r denote a random number from the uniform distribution on the 
interval 0 to 1. If r = F(x), then x = F(r) where F(r) ^ denotes 
the inverse transformation of F(x). If X^ is exponentially distributed 

The transformation of a random variable X into a uniformly distributed 
random variable r = F(x) is a well known result in probability theory. 
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with parameter X , then 



F(x.) = 1-e or if F(x.) = r., 

then r. = 1-e ^^i; 

1 

"Xx 

thus e i = 1 - r but if r. is a realization from the uniform dis- 

1 1 

tribution on the interval 0 to 1, then 1 - r. = r.’ is also a realization 
on the same interval. Thus this may be rewritten 



-Xx. , 

e 1 - r . 

1 



solving now for x^, 



k 

If X = S X. , 
i=l ^ 

X ; thus 



x. = - 1/X In r . * . 

1 1 

then X has an Erlang distribution with parameters k and 



k 

X = - 1/X E In r. * ; which may be rewritten 
i=l ^ 



X = - lA In ( TT r. 

1=1 ^ 

F(r)‘^ = - lA In Q 




where 



k 

TT r 
i=l 



.0 



with r 




3 . Test for Goodness of Fit by Chi-square Criteria 

Let x^ be. a realization from the sampled population; let Oj denote 
the number of realizations that fall into class j for j = l,2,.,.,n, 
where n is equal to some arbitrary number of subdivisions of the range 
of the assumed distribution of X^ ; let Ej denote the expected number 
of realizations in class j under the assumed distribution; then the 
stat istic 



n (E. 



X = Z 
i=l 






E. 

J 



has an approximate chi-squared distribution with parameter (n-m-l)/2, 
where m denotes the number of assumed distribution parameters estimated 
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m 




with sample statistics, n-m-1 is called the degree of freedom associ- 
ated with X; for a choosen significance level the value of X can be 
compared with the appropriate table value such as found in Appendix ^ 
of [5]; if X is greater than the table value, then the null hypotheses 
that the came from the assumed distribution is rejected; it is 
accepted otherwise . 

4. Estimation of Gamma Parameters by the Method of Moments 

Let T be a random variable from a gamma distribution with parameters 

k and X. Let denote the expected value of T, T denote the T sample 

2 2 

mean value, Orj, denote the variance of T, and denote the T sample 

"" - 2 

variance as calculated by l/(n-l) E (T.-T) , where n is the number of 

i=l ^ 

T sample values. Since 



= kA, and 

2 - ^a2 

= k/X , 



A A 

it can be shown that k and X, estimates of k and X, 
be calculated by 



A 

k 



T^/S^ and 






T/S^. 



respectively, can 



if the mean \ji and variance are respectively estimated by the sample 
mean and variance. 

5 . Estimation of Gamma Parameters by the Peak Value Method 

Let T be a random variable from a gamma distribution with parameters 
k and X. Let denote the expected value of T, T denote the T sample 
mean value, p denote the maximum value, i.e., the peak, of the relevant 
gamma density function, and m denote the mode of the distribution of T, 

A /y 

i.e., the T value at which p occurs. Let |_i^, p, k* and X* denote 
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estimates of |j^, p, k and X, respectively. Since p occurs at m and 
since m = (k-l )/\5 p can be evaluated as 



,,k-l 1-k 

„ = (k-l) e 

P r(k) 

Now by applying an approximation formula for T (k) as found in [6], 

i.e. , 

j — k-1/2 

F(k) = V2n e k for k large and not an integer, 

p can be approximated as 



p^— 

v2rrk 

k-1 -1 

Since for large k, (1-1/k) can be approximated by e , p can be 
approximated as 

( 1 ) ■ 

/s/2nk 

A A 

To estimate k and X, and p will be used. Since = k/X and since 
p can be approximated by formula (I), by combining these equations and 
solving for k, it is easily seen that 

k ^ 2TTp^^^. 

Now, knowing an approximation for k, an approximation of X can be 

obtained from X = k/|a^. Now, by letting |j^ == T and p = max Oj/nw, 

where max 0. is the number of T realizations in the modal observation 
J 

class of the T sample frequency histogram, n is the number of T 
realizations in the sample, and w is the width of the modal class in 
the frequency histogram, k* and X* can then be evaluated as 

k* = 2np^T^ 

\* = k*/T. 
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